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Abstract 



We outline reasoning that led to the current theory of quasars and look 
at George Contopoulos's place in the long history of the N-body problem. 
Following Newton we find new exactly soluble N-body problems with multibody 
£^ \ forces and give a strange eternally pulsating system that in its other degrees of 

£SJ ■ freedom reaches statistical equilibrium. 



OO 

Introduction 

^ ■ 

Q-f This meeting is held to honour George Contopoulos for his great contributions to 

O ' dynamical systems theory and the N-body problem. I shall pay my tribute to him 

£3 , in three parts, 

CO 



I Placing his contributions in the proud history of those who have made major 
contributions to the N-body problem. 



II Since nothing but the best is good enough to honour George I present to him 
a copy of my best paper (Lynden-Bell, 1969) and include here a resume of its 
arguments that led to the current theory of quasars. There are questions my 
paper raised 29 years ago which are still unexplored. 

Ill With Prof. Ruth Lynden-Bell (my wife) I present our new extraordinary 
N-body problems which we solve for all initial conditions. These problems 
can also be solved in quantum mechanics when the hyper-keplerian potential 
energy is 

V = -NZe 2 /r , (1) 



where 



JV 
=1 



M = Y,mi (3) 
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i.e., r is the mass- weighted-root-mean-square radius of the N-body system. I 
quote here the energy and degeneracy of the nth quantum state but we shall 
publish derivations elsewhere. 

2Mh~ 2 (NZe 2 ) 2 

E n = ^ V (4) 

[2n + 3(iV-2)] 2 
the degeneracy of this N particle state is 

for N = 2, g reduces to n 2 and the energy reduces that of the hydrogen atom 
for which M = m p + m e . To see this, Z is replaced by \{m v m e ) 1 l' 1 jM when r 
is replaced by the separation of the electron from the proton. We have written 
the coefficient of the potential energy in the clumsy form NZe 2 so that the 
analogy to hydrogenic atoms can be readily seen by physicists. 

1 Contributions to the N-body problem (excluding Agatha 
Christie's) 

The N-body problem probably started with Newton although Hooke would un- 
doubtedly dispute it since he seems to have conceived the idea independently but 
had not the mathematical ability to work out its consequences. As Chandrasekhar 
(1996) has shown, Newton's Principia (Newton 1687, Cajori 1934) has much to 
teach us even today (see Lynden-Bell & Nouri-Zonoz, 1998). Recent studies of the 
Portsmouth papers have shown that Newton developed most of the perturbation 
theory that was hitherto attributed to the mathematical astronomers of the 18th 
and 19th centuries. Newton's method was to store up the momentum generated by 
perturbations and then deliver it as an impulse that changed the motion from one 
ellipse to another. This of course gives him the equations for the variations of the 
orbital elements which are the meat of perturbation theory. My brief resume of the 
N-body problem's history is: 

Newton 1687 Orbit Theory and the general solution of 

the first extraordinary N-body problem 

Laplace 1795 Perturbation Theory for near circular 

orbits 

Poincare 1892-99 Topological Methods 

Whittaker 1913 (1959) Adelphic Integrals as series 

Contopoulos 1956- Third Integrals and Chaos 

Kolmagorov-Arnold-Moser Invariant Tori & Arnold diffusion 
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To these theoretical studies we must add the numerical computation of the N- 
body problem and here Aarseth's name stands out as a persistent pioneer exploring 
this problem (Aarseth, 1974) although many others have contributed, especially 
Heggie (1975) through his work on triple interactions. Both in globular cluster 
theory and in dynamical systems Henon's work stands out for its beauty (Henon, 
1961, 1969, 1974) while Antonov (1962) was responsible for a fundamental advance 
in the understanding of gravitational thermodynamics later popularised and ex- 
tended to negative specific heats and gravitational phase transitions by the author 
(Lynden-Bell & Wood, 1968; Lynden-Bell & Lynden-Bell, 1977)and by Thirring 
(1972). Betteweiser & Sugimoto (1984) were responsible for giving the gravother- 
mal instability a delightful new twist in their discovery of the inverse gravothermal 
catastrophe that leads to giant thermal oscillations. But let me return to what 
George Contopoulos taught me at our many contacts since 1961. 

To set the scene I had written my thesis in 1960 which contained a new derivation 
of what potentials had local first integrals of the motion besides the energy and the 
angular momentum about the axis. The main part of the work was the derivation of 
these different classes of potential while other parts of the thesis contained the time 
dependent evolution of accretion disksQ and a first attempt to apply Jeans's (1928) 
gravitational instability to make a theory of the spiral structure of galaxies. The 
beliefs of those times are well illustrated by the first edition of Landau & Lifshitz's 
book on classical mechanics; either a dynamical system was separable and integrable 
or it was ergodic (by which was meant that almost all orbits visited all volumes of the 
phase space accessible under the energy constraint). Having classified the special 
forms of potential that had local integrals I expected that most other potentials 
would show ergodic behaviour. From the inequality of the z and R dispersions 
of the stars in the Galaxy it was clear that there must be another integral other 
than E and h for the Milky Way so I had begun trying to fit Eddington (1915) 
(now called Stakle) potentials to galaxies. It was quite shattering when at the 1961 
IAU general assembly in Berkeley, George Contopoulos (1960) showed that orbits 
in most smooth potentials behaved as though there were third integrals. Suddenly 
the special interest of the special potentials fell away — they were not the only 
systems with 3rd integrals, merely those for which we knew the exact analytical 
form of those integrals. They seemed now to be mathematical curiosities rather 
than systems fundamental to the dynamics of real galaxies. 

Three years later George organised a very instructive IAU symposium (No. 25) at 
Thessaloniki on the Theory of Orbits in the Solar System and Stellar Systems. Here 
he brought into contact the celestial mechanics fraternity, with their long history of 
analytically calculating orbits in the solar system by perturbation theory, with us 
new boys who were attempting to understand the statistics of the orbits in the more 
complicated potentials of galaxies; George Contopoulos (1965, 1966) here taught us 
that many of the problems were common to both fields and showed how fertile it 
was to bring different communities who knew different things to the same conference 
- his wide interests have made him especially good at that throughout his life and 

1 they got that name only later 
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this 1998 conference is no exception. 

For brevity, I shall skip contacts at Besangon on the N-body problem where 
George presented Poisson Bracket series for third integrals and we were introduced 
to Lie series. 

In 1973 at Saas Fee George gave lectures in which he introduced me to the 
wonders of modern dynamical theory - topological methods incomplete chaos and 
the KAM theorem. It opened my eyes to so much that was new to me that I retreated 
back to more directly astronomical topics preferring the contact with astronomy to 
the unchartered seas revealed by this new alliance between computers and topology. 

Two years ago at Salsjobaden in a conference on the Dynamics of Barred Spirals, 
George again broke open a new field (Contopoulos, 1997). His invariant dynamical 
spectra (described also in his contribution here) taught us how to measure and 
classify chaos, even complete chaos! 

I have picked out a tiny fraction of George Contopoulos's work (1975) and men- 
tioned things I learned from our direct contacts. He will no doubt deduce that I 
am not a very attentive pupil but it would be mean not to mention a lovely paper 
on the light distributions of elliptical galaxies (Contopoulos, 1956) because it is a 
beautiful work to which I constantly have to refer my astronomical colleagues! 

The essence of this paper can be deduced by the following argument. Consider a 
spherical galaxy with any radial light profile. Now flatten its density distribution by 
linear contraction along any axis. This contraction can be resolved into one along 
and one perpendicular to the line of sight. The one along makes no difference while 
the one perpendicular flattens the circular distribution of observed light into one 
stratified on similar ellipses. If a further contraction is made along another axis we 
can apply the same argument again since ellipses contracted along any direction 
remain ellipses. So we arrive at George's beautiful theorem that if the density 
distribution of an elliptical galaxy is stratified on similar concentric ellipsoids then 
the light seen will be stratified on similar concentric ellipses whatever the orientation 
of the galaxy to the line of sight. 

2 Background to the Accretion Disc Theory of Quasars 

My own best work is "Galactic Nuclei as Collapsed Old Quasars" written in 1969. 
Then the discovery of quasars by Schmidt using Hazard's accurate position for one of 
Ryle's radio sources was still recent and quasars themselves were enigmatic objects 
more especially so because even the brightest 3C273 and 3C48 too did not seem to 
be associated with clusters of galaxies. 

No-one then knew that Michell in his wonderfully percipient paper of 1784 had 
predicted both giant black holes and how they would be discovered! Even the 
name black hole only came into general use in 1970! In my 1969 paper I refer 
to "Schwarzschild throats". Laplace's translation of Michell's work into French 
(without attribution!) was not common reading among astronomers either. 

Among the modern works on quasars as accretion discs, priority goes to Salpeter's 
fine 1964 letter to the Astrophysical Journal. Turning against the then common view 
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that quasars were not associated with clusters of galaxies he worked out the conse- 
quences of a large black hole moving through a galaxy and accreting according to 
the Hoyle & Lyttleton formula. He derived the power emitted per unit accretion 
rate by considering the binding energy of the last stable circular orbit and deduced 
a number of consequences of such black holes accreting as they wandered through 
the interstellar gas of a galactic disc. Five years elapsed before I wrote my paper. 
Originally unaware of Salpeter's (1964) note I luckily learned of it before the proofs 
came and so was able to add a sentence and a reference to his work. My aim was 
to show that the very small nuclei already known in the centres of galaxies were 
likely to be stars gathered around the giant-black-hole remnants of quasars. At the 
time, 1969, we already knew that the Optical Violently Variable or OVV quasars 
could change by a magnitude from one night to the next. Geoffrey Burbridge (1958) 
had been insistent that the giant radio sources needed 10 61 ergs in fast electrons 
and magnetic field, while Ryle (1968) had emphasised that quasars would not be 
distinguished from such sources by radio measurements. 

Now 10 61 ergs weigh ^10 7 M Q . If one entertained the idea that these ergs came 
from nuclear energy then the 1% mass conversion efficiency of nuclear burning means 
that 1O 9 M are needed. However putting 10 9 M Q within the light-variation-time 
length-scale of 10 light hours gives a gravitational binding energy of 10 62 ergs - on 
such a hypothesis 10 62 ergs of gravitational energy would have been lost, all in order 
to burn 10 9 Mq of hydrogen into Helium and thereby get the mere 10 61 ergs needed. 
This shows that in assuming nuclear power we nevertheless conclude that most of 
the energy comes from gravity. So the nuclear idea is not sensible and we should 
assume a preponderant gravity power and a somewhat smaller mass ~ 10 8 Mq. If 
conversion of mass into radiation is not 100% efficient quasars must leave behind 
massive remnants of > 10 7 — 1O 8 M and because they have radiated their binding 
energy they have insufficient energy to re-expand. Since the masses are far beyond 
the Chandrasekhar limit there are no other final resting places other than giant 
black holes. Turning to the numbers of quasars derived by Sandage and estimating 
possible lifetimes, I deduced 

Number of clusters Number of dead Number of 

of galaxies < quasars < galaxies 

Thus the nearest dead quasar must be nearer than M87 and there may be as many 
dead quasars as there are massive galaxies. How could we hide dead quasars of 
1O 8 M when they still gravitate? They would naturally be centres of attraction for 
stars so it is natural to find such a body at the centre of an exceptionally dense 
region. Galactic Nuclei then became the obvious candidates so I looked at the 
Galaxy, the Magellanic Clouds, M31, M32, M81, M82, NGC4151, M87, etc. and 
estimated possible black hole masses from the 1969 data on their nuclei, many of 
which were due to pioneering work by Merle Walker, see Figure 1. I also drew on the 
accretion discs of my thesis and, finding the gaseous viscosity too low, I estimated 
a magnetic viscosity. This was based on the shearing of the disc causing magnetic 
reconnection and continual flaring above the disc. Indeed I found that the protons 
got most of the energy as they more readily achieved "runaway" . Particle energies 
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Figure 1: Comparison of my 1969 estimates of black hole masses O with modern 
determinations * and less precise modern estimates •. Evidently the early estimates 
were somewhat over enthusiastic! 



up to 10 eV were readily generated and hard emission would result when this hit 
the disc material. While the energy was primarily dissipated into such fast cosmic 
rays they would collide with the disc and heat it to temperatures T oc r~ 3 / 4 for 
r 3> 2GMo/c 2 . Adding together such black body rings of emission, I got the disc 
spectrum S u ~ v 1 ^ exp — (hv/kT max ) where T max the maximum temperature in 
Kelvin was 6.6 x 10 4 F^M 7 here F_3 is the mass flux in units of lO _3 M /yr 
with Mj the black hole's mass in units of 10 7 M Q . I did not estimate how much 
hard emission would come from the initial collisions of the cosmic rays with the disk 
but a 10 13 eV cosmic ray is certainly capable of emitting hard 7 rays at its first few 
collisions. Even today, 29 years later, I think this model deserves more attention 
as a serious rival to the currently popular advection models. 

The following year Jim Bardeen wrote a particularly fine paper which showed 
how accretion would spin up a Schwarzschild hole and after a finite mass was ac- 
creted leave it growing as a near-limiting Kerr hole of significantly greater efficiency. 
I gave a paper on these models at the 1970 Vatican Symposium on the nuclei of 
galaxies and in 1971 re-estimated the luminosity functions of quasars and mini 
quasars by developing the C~ method. That year Ekers was stimulated to look 
with higher radio resolution at the Galactic Centre (Ekers & Lynden-Bell, 1971) 
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and I reviewed the then known data on it with Rees (Lynden-Bell & Rees, 1971). 
A year or two later attention turned to lower mass black holes with the discovery 
of many X-ray binaries by the UHURU satellite. Papers by Pringle & Rees (1972) 
and Shakura & Sunyaev (1973, 1976) applied such ideas on a smaller scale and with 
Pringle I applied them (1974) to star formation both with and without magnetic 
fields. In 1978 I introduced the thick Kerr and Schwarzschild vortices in the hope 
of getting a more natural collimation mechanism than that of Blandford & Rees 
(1974) but the very narrow jets are still inadequately understood. 

3 General Exact Solution to an Extraordinary N-body 
Problem 

I now return to the N-body problem and the little known fact that in Principia 
Newton (1687) solved an N-body problem in which every body attracts every other 
one and he solved it for all initial conditions! 

His was the first of the class of extraordinary N-body problems which Ruth 
Lynden-Bell and I have been studying. Newton took the force between two bodies 
i and j to be F = kmirrij(xj — Xj). To get the total force in particle i he summed 
over j and since the j = i term is zero we may sum over all j to obtain 

Fi = J2 F ij = kmiM(x - Xi) (6) 

j 

where x is the position vector of the centre of mass which of course moves uniformly 
in a straight line and M is the total mass of the system. Thus with this linear mass- 
weighted law, that Newton would never have ascribed to Hooke, the total force on 
the i th body is directed to the centre of mass and proportional to the distance from 
it. Therefore Newton found that each body describes a centred ellipse about the 
centre of mass which itself moves uniformly. This completes Newton's solution. In 
his case the potential energy is 

V = -\K ^mimjfc -xj) 2 = -|*M^ m;(x; - x) 2 = -\kM 2 r 2 . (7) 

i < j i 

Generalising some work on statistical mechanics by Ruth Lynden-Bell we were 
led to consider the dynamics of N-body systems with the more general potential 
energy V = V(r) where r is given above, (cf. equation (2)). We define a mass 
weighted radius r in 3iV dimensions by 

r= WM (xi " x) ' Vm (X2_x) ' vir ( xjv ~ x v ' (8) 

so the first 3 of the N coordinates tell us where particle 1 is, the next 3 where 
particle 2 is, etc. Notice that |r| is the r we defined previously. Equations of motion 
of the particles in centre of mass coordinates then lead directly to the equation 

Mr = -V'{t)t (9) 
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where r = r/r is the unit radial vector in 3N space. 
One readily sees that 

r a rp - rpr a = , (10) 

so 

r a r/3 - rpf a = L a/3 = —Lp a = const . (11) 

Furthermore 

L 2 = \L a/3 L a p = \ (r a ff3 - rpr a ) (r a r^ - rpr a ) = r 2 (r) 2 - (r • r) 2 (12) 

where L 2 is the constant defined by the first equality. 

The energy in centre of mass coordinates is given therefore by 

±Mr 2 + V(r) = E = \M(r 2 + L 2 r~ 3 ) + V(r) . (13) 

This determines r{t) as a periodic function if E < so there is no violent relaxation 
in these systems and they vibrate eternally. 
Differentiating (13) we find 

M(r - L 2 r" 3 ) = -V'{r) . (14) 

This is the same equation of motion as that for the central distance to an object in 
planar motion which angular momentum L about a centre of force with potential 
V{r). It is natural to imagine such a planar orbit and to invent an angle <f> such 
that 4> = at some pericentre and 

r 2 4> = L, (15) 

we may then imagine an orbit in two dimensional polar coordinates r, <p and following 
Newton we shall cling to the geometry by eliminating the time in favour of <p. Now 

.. d 2 . A , d ( A Ldr\ A L 2 d 2 i , . 

r = df 2{rr) = dA rY+ -^) =rY+ ^w (16) 

where two terms in rLr~ 2 dr /dtp cancel at the last step. Inserting this result into our 
equation of motion (9) and using (14), we deduce the wonderfully simple equation 

d 2 f/# 2 + r = , (17) 

whose solution is 

r = Acos</> + Bsin0 (18) 

where A and B are constant 3N- vectors which obey A 2 = B 2 = 1 and A • B = in 
order that r should be a unit vector for all (j). Three further constraints on A and 
B follow from the fixed centre of mass. They are detailed in our paper but need 
not concern us here. 

We now have the general solution, the centre of mass moves uniformly in a line 
and the particles pursue orbits about it of the form 

r = r((p)(A cos <j> + B sin <f>) , (19) 
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where r(</>) is the form of the two dimensional orbit governed by equations (13) and 
(15). These can be integrated explicitly for the Isochrone potential V oc k/(b + s), 
s 2 = r 2 + b 2 and for the Kepler and harmonic oscillator potentials. For the Kepler 
case r((f>) = £/(l + ecos^) so the solution is of the pleasing form 

r = £(l + ecos0) -1 (Acos</> + Bsin0) . 

If we concentrate on the particle i, we find its orbit lies in the plane perpendicular to 
Aj x Bj where i denotes the three components corresponding to particle i. Taking 
x, y coordinates in that plane and eliminating <p we find that the orbit is quadratic. 
If e were zero it would be a central ellipse, while if |Aj| and |Bj| are equal and 
orthogonal it gives a Keplerian eccentric ellipse. In the general bound case the 
ellipse has neither its centre nor its focus at the centre of mass r = 0. These 
systems obey the equilibrium Virial theorem in the form 2T — rV = 0, so for the 
hyper-keplerian case V oc r _1 it takes the more familiar form 2T + V = 0. 

One may work out the microcanonical statistical mechanics and find that 
E = -§(iV - 2)kT so that the heat capacity C = — |(iV - 2)k which is clearly 
negative as for other gravitating systems (Lynden-Bell & Wood, 1968) and black 
holes. If V takes the form 



V 



oo r < b 

-kM 2 /r b <r < R 

oo r > R 



corresponding to a gravitating system which cannot get too small or too big then 
a Canonical ensemble is possible and the negative specific heat region of the mi- 
crocanonical ensemble is replaced by a giant first order phase transition as in our 
earlier model (Lynden-Bell & Lynden-Bell, 1977). 



3b Generalisation 

We may extend these extraordinary N-body problems by taking V to be of the more 
general form 

V = V {r) + r- 2 V 2 {r) 

the only restriction on the second term being that it scales under expansion as r~ 2 . 
Those familiar with separable systems in 3 dimensions will know that for such poten- 
tials (imh 2 + V2) is constant along an orbit where for that case 
V2 = V2(9,4>) and h = r x v. The generalisation to 3iV dimensions is the first 
integral \ML 2 + V^(r) = \MC? say (note that due to the V2 term, C 2 does not 
have to be positive). 

The energy equation now reads 

E = \Mi 2 + V Q + r~ 2 V 2 = \M{f 2 + L 2 r~ 2 ) + V + r~ 2 V 2 = \M{r 2 + C 2 r~ 2 ) + V 

so the r motion pulsates for ever as before. These systems show no violent relaxation 
in their breathing mode which pulsates (or evolves E > 0) independently of the 
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complication of the r motion. Since V2(r) is still free to choose, that motion can be 
as complicated as we like to make it. Defining a new time r by d/dr = l/r 2 d/dt 
the equations of motion for r as a function of r are totally independent of the r 
motion, having a reduced Lagrangian system of their own in r-time. An interesting 
case to consider is the statistical mechanics of a "hard cone" gas in which V2 is 
large and repulsive only in very small regions where two particles are nearly in the 
same direction as seen from the mass centre. This corresponds to the small hard 
sphere gas so beloved of textbooks. Carrying out that statistical mechanics, which 
is totally independent of any r motion that may be going on, we obtain a new 
system at equilibrium in its r coordinates but pulsating or evolving in r. 

We have shown (Lynden-Bell & Lynden-Bell, 1998) this equilibrium to be best 
described in terms of the peculiar velocity Vj relative to a "Hubble flow" i!f(xj — x) 
where H = r/r that is 

Vj = ij - x - H(xi - x) 
/(Vj,x - x) oc exp 

Thus the distribution is Maxwell-Boltzmann relative to be mean Hubble flow with 
a temperature proportional to r~ 2 (t) and the profile is gaussian with a dispersion 
proportional to r(t). It is notable that the 'equilibrium' of the r coordinates is main- 
tained throughout the pulsation just as the Planck distribution of cosmic black-body 
radiation in the Universe is maintained without interaction during the expansion of 
the Universe. Thus whether the relaxation to equilibrium of the angular coordi- 
nates is longer than or shorter than the pulsation time of r is not relevant because 
'equilibrium' once attained is maintained throughout the pulsation, it does not have 
to be recreated as each radius r is attained. 
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